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ABSTRACT 

A new field of numerical astrophysics is introduced which addresses the solution of 
large, multidimensional structural or slowly-evolving problems (rotating stars, interact- 
ing binaries, thick advective accretion disks, four dimensional spacetimes, etc.), as well 
as rapidly-evolving systems. The technique employed is the Finite Element Method 
(FEM), which has been used to solve engineering structural problems for more than 
three decades. The approach developed herein has the following key features: 

1. The computational mesh can extend into the time dimension, as well as space — 
generally only a few cells deep for most (fiat-space) astrophysical problems, but 
throughout spacetime for solving Einstein's field equations. 

2. When time is treated as a mesh dimension, virtually all equations describing the 
astrophysics of continuous media, including the field equations, can be written in 
a compact form similar to that routinely solved by most engineering finite element 
codes (albeit for nonlinear equations in a four-dimensional spacetime instead of 
linear ones in two or three space dimensions): the divergence of a generalized stress 
tensor equals a generalized body force vector, both of which are functions only of 
position, the state variables and their gradients. 

3. The transformations that occur naturally in the four-dimensional FEM possess both 
coordinate and boost features, such that 

(a) although the computational mesh may have a complex, non-analytic, curvilinear 
structure, and may be adapted to the geometry of the problem, the physical 
equations still can be written in a simple coordinate system that is independent 
of the mesh structure. 

(b) if the mesh has a complex fiow velocity with respect to coordinate space, the 
transformations will form the proper advective derivatives, automatically con- 
verting the equations to arbitrary Lagrangian-Eulerian. 

4. Only relatively simple differential equations need to be encoded. The complex dif- 
ference equations on the arbitrary curvilinear grid are generated automatically by 
the FEM integrals. A different integration method must be used for equations of 
odd and even order. 



This first paper concentrates on developing a robust and widely-applicable set of tech- 
niques using the nonlinear FEM and presents some examples. The second paper in this 
series will deal with making the method fast and efficient, so that large, astrophysically- 
interesting computational meshes can be employed. 

Subject headings: methods, numerical — hydrodynamics — magnetohydrodynamics: 
MHD — relativity — stars: rotation 
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1. Introduction 

The first problem to be solved with the techniques 
of numerical astrophysics was the structure and evo- 
lution of stars — an "implicit" proble m that involves 
a sta t ic or slowly-evolving struc ture ( Chandrasekhar 



1957 ; AUer fc McLaughlin 1965 ). Its solution consists 



of determining the state variables of the fluid (den- 
sity, temperature, pressure, flux of radiation, compo- 
sition) at each radius in the stellar interior, and is 
obtained by relaxing a large set of coupled nonlin- 
ear difference equations (derived from the differen- 
tial equations of stellar structure) along with bound- 
ary conditions at the stellar center and surface. Ra- 
dial stellar structure is only a one-dimensional prob- 
lem and, while once considered difficult and CPU- 
intensive, now is solved easily on personal comput- 
ers (PCs). Since then, many other fields of numeri- 
cal astrophysics have been developed: "explicit" hy- 
drodynamic simulations of explosive and jet phenom- 



ena (Norman 1997); N-body and smooth particle hy- 
drodynamics (SPH) of discrete or semi-discrete sys- 



tems of particles (Dubinski & Hernquist 1997; Mon 



ighan 1992); Monte Carlo simulations of radiation 
flow ([Leahy 1997|; [Park fc Hong 1998|); etc. AU have 



matured to the point where the solution of three- 
dimensional, time-dependent problems is not uncom- 
mon. 

Ironically, however, only modest progress has been 
made in extending the original implicit problems into 



sev eral dimensions : — rapidly rotating stars, evolving 



an d interacting binaries, detailed accretion disk struc- 
ture and evolution, etc. There are several important 
reasons for this. Firstly, the geometry of these sys- 
tems is unknown until the problem is solved. For 
example, the shape of the outer surface of a rotating 
sta: 



or (possibly thick and advective) accretion disk 
wil l be part of the solution, and the shape of a rapidly 



rotating stellar core may have a different oblateness 
(or even prolateness) from that of the outer envelope. 
No numerical method capable of operating in nearly- 
arbitrary geometries has been applied extensively in 
astrophysics. Instead, one either has assumed spher- 
ical symmetry and treated only slowly-rotating, per- 
turbation problems (Kippenhahn & Thomas 1970), 
or has assumed that the isosurfaces of the state vari- 
ables are coincident (which implies rotation on cylin- 
ders via von Zeipel's theorem) and again solved an es- 
sentially one-dimensional, or limited two-dimen sional 
problem ([Eriguchi fc Miiller 199"ll ; [Clement 199^ ). Re- 



cently, some progress for two- and three-dimensional 
stellar models has been made using a multi-domain 
approach and treating the stellar surface as a discon- 
tinuity ( Bonazzola, Gourgoulhon, fc Marck 1998[ ). 

Secondly, even if a general geometrical method 
applicable to large numbers of two-dimensional and 
three-dimensional problems could be developed, the 
current stellar structure methods for solving the im- 
mense system of simultaneous nonlinear equations 
would take a prohibitively long amount of CPU time 
and memory. For example, a relatively modest prob- 
lem with 256'^ grid points and ten state variables at 
each point would generate a banded matrix 10*'^ x 
10®-^ in size, taking up at least 10^'*-^ bytes (0.9 
PB) for the non-zero elements. Direct inversion tech- 
niques, similar to the Henyey method commonly used 
in stellar structure (which take the handedness into 
account), would take about a thousand years to invert 
this matrix once on a large parallel supercomputer 
like the Cray T3D, with perhaps lO'' or more such 
inversions necessary for a complete stellar evolution 
modelQ 

Fortunately, there exist techniques for solving both 
of these problems that are well developed and have 
been in use in the engineering field for many years 
(although it is still rare to see both used at the same 
time). The Finite Element Method (FEM), intro- 
duced more than four decades ago and the preferred 
method of treating multidimensional structural en- 
gineering problems since the late 1960s (Zienkiewicz 
1977), approximates objects as distorted lattices of 
small structural members called elements. For solv- 
ing large systems of coupled equations generated by 
such grid problems, the multigrid method was intro- 



duced in the 1970s (Brandt 1977) and is now begin- 
ning to be used in astrophysics as well (Truelove ei 
1998|; Norman 1998). This approach dispenses 



al. 



with the large matrix, cleverly reaching a solution af- 
ter a few sweeps of the mesh. Together, these tech- 
niques promise to make multidimensional astrophysi- 
cal structural problems possible, and bring the time to 
solve them within an order of magnitude or so of that 
for explicit problems. When coupled with the contin- 
ued expected increase in speed of computers over the 
next few decades (which has averaged about a factor 
of 2 every 2 years for the past 20 years), it is not 



^For this calculation it is assumed that the storage is propor- 
tional to the bandwidth B (~ 10^'*) of the matrix times its 
length L and the time to invert is proportional to L. 
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inconceivable that three-dimensional structure prob- 
lems will soon be solved routinely on future models of 
PCs and that four-dimensional problems will become 
commonplace on supercomputers. 

This first paper deals with the development of 
the general geometrical method for solving multi- 
dimensional structure and evolution problems. At 
this stage the speed and efficiency of execution of 
the method will not be a concern; the focus will be 
only on producing a robust and widely-applicable set 
of useful techniques. Our goal will be to determine 
the essential features of most astrophysical systems of 
equations and the geometrical demands they place on 
the numerical method. These properties will then be 
encoded at the outset, ensuring some measure of gen- 
erality. The second section describes the set of equa- 
tions that can be addressed with nonlinear astrophys- 
ical finite element analysis and develops a method for 
solution on fixed (non- moving) grids. In section 3, 
this then is generalized to include situations where 
the positions of the grid points are part of the solu- 
tion and the grid can change with time. Finally, tests 
and examples are given, using the author's code, in- 
cluding rotating polytropic star models. 

2. The Basic Four-Dimensional, Nonlinear Fi- 
nite Element Method on Fixed Grids 

This section describes the techniques used in the 
author's computer code, entitled GENRAL, for solv- 
ing general astrophysical problems. It utilizes the 
techniques of finite element analysis (FEA) — in use 
in the field of engineering for some time — but gener- 
alizes them to nonlinear equations in four dimensions, 
instead of linear equations in two or three dimensions. 
For this initial development, it is assumed that the 
coordinates of the grid points (or "nodes") at which 
the variables are evaluated do not change while the 
solution is being computed. 

2.1. General Form for Equations of Contin- 
uum Astrophysics 

Appendix A shows that the differential equations 
of continuum astrophysics in curved spacetime can be 
cast into the generic form 



= (T^^iVw, w, x)^g)^p 
— FqCVw, w, x) \J—g 



where [to]" = is a generalized solution vector hold- 
ing all of the w = 1, y unknown state variables; 5Jg 
is a generafized residual for each of the q = l,...,y 
equations (which will be forced to zero through nu- 
merical relaxation techniques); Tq and Fq are, respec- 
tively, the generalized stress tensor and force vector 
for these equations; and g is the determinant of the 
metric tensor g ([ — h -I- -signature) In a similar 
manner, the boundary conditions on these equations 
can be cast as 

3?, = w, x)^g\^,S^ii 

-f • Vio, w, a;) 

-./r(Vl(?, W, X)^f^ 

= (2) 

for each of the r = 1, boundary conditions. In 
general, R^V since, depending on the highest order 
derivative in 3?^, there may be 0, 1, or 2 associated 
boundary conditions. The last term in (H) is adequate 
for handling Dirichlet, Neumann, and mixed bound- 
ary conditions, such as the radiative condition at a 
stellar or accretion disk surface. The first two terms 
are necessary for including constraints on field equa- 
tions. S is the projection tensor along the boundary 
9r2 and orthogonal to the boundary normal n 



S = n ® n 



(3) 



where n • n — —1, n • S ~ 0, and S • S ~ S, "(8)" 
is the outer (dyadic) product, and "•" is the inner 
(scalar) product. S causes the divergence in equation 
(H) to be performed on the boundary only and the 
derivative normal to the boundary to be only first 



Thrniiffhniit ih, 



paper the notation of Misner, Thorne, & 



Wheeler (1973) is used, with Greek letters indicating coor- 
dinate indices in four-dimensional spacetime (« = 0, 1, 2, 3) 
and Latin "integer" letters indicating three-space indices only 
(i = 1,2,3). The comma denotes ordinary differentiation with 
respect to the coordinates, while a semicolon will denote co- 
variant differentiation. Repeated indices indicate summation 
over the entire range of those indices (the Einstein summation 
convention), so that 



9 9iil3,j 



dx-i 



(1) 



A raised index indicates contravariant properties of the tensor 
and a lowered index indicates covariant properties. Note that 
w'" is written as a contravariant vector with a raised index, like 
the coordinates x"; this is partly for convenience (to facilitate 
the summation convention) and partly to draw attention to the 
variables as generalized coordinates of the system. 
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order in d/dx^. An alternative form for (Q) is 

= [if (Vii;, w, x)^]^^S^'p 

^[h";{S-Vw, w, 

-/V(Vu;, w, x)y/^ 
= 

which also has no second derivatives in d/dx^ . 

It is important to note that equation ([T]) includes 
not only structural and steady problems, but also 
evolving ones as well. For these cases, in addition to 
having three spatial coordinates, the computational 
grid can extend into the fourth (time) dimension, pos- 
sibly from the initial time step or hypersurface to the 
final one. While certainly increasing the computa- 
tional and memory load on the computer, this ap- 
proach will have distinct advantages over conventional 
approaches to initial- value problems. 

The boundary conditions are not necessarily com- 
pletely described by equation (||). While it is well 
known that, physically, the boundary of a boundary 
is zero {ddVl — 0), computationally one often intro- 
duces sub-boundaries by truncating the mesh or im- 
posing symmetries on the problem. These conditions 
produce right-angle kinks in the boundary, where n 
suddenly rotates by 90° and the boundary conditions 
abruptly change. Such boundary corners occur, e.g., 
where the t = Q initial hypersurface intersects the 
world line of a stellar surface or (in the case of ax- 
isymmetry or plane symmetry) where the stellar sur- 
face intersects the symmetry axis or plane. For exam- 
ple, when solving Maxwell's evolutionary equations 
on the four-dimensional domain f2 under such con- 
ditions, they will be bounded on the t = Q portion 
of dVt by the initial value (solenoidal and Coulomb) 
constraints; these will be bounded further at an exter- 
nal stellar 2-surface dd^; and these may be bounded 
still further by the rotation axis or equatorial plane at 
edges dddfl. Therefore, additional equations, similar 
to (H), with successive projection of the first two terms 
into the sub-boundaries of lower dimension, may be 
needed until one reaches the zero-dimensional ddddfl 
(endpoints of line segments) where the conditions be- 
come simply 

= -!r{Vw, W, X) = Q (4) 

In the examples in this paper all boundary conditions 
are of the simple form (Q) , but in general astrophysical 
situations the form (El) will be needed. 



2.2. Continuous Solution: The Element Mesh 

Formally, in the finite element method (FEM) the 
computational domain SI is subdivided not into nodes, 
but into sub-domains {5il) called "elements" — simi- 
lar to "zones" or "cells" in the finite difference method 
(FDM). (Nodes come later, and then only to facilitate 
the element process.) The elements are constructed 
in such a way that each function is continuous over 
the entire domain, but its derivatives are only piece- 
wise continuous; i.e., is continuously differentiable 
only within each element. The cc" are treated in the 
same manner; they also are continuous across element 
boundaries with no spatial "gaps" between elements. 

Although a variety of generic element shapes can 
be used, the most common are triangular and quad- 
rangular. Of course, these assume higher-order shapes 
in three and four dimensions (i.e., equilateral trian- 
gles, tetrahedra, simplices [hyper-tetrahedra] ; squares, 
cubes, and hypercubes), but they still shall be re- 
ferred to here as the triangular and quadrangular 
classes. In GENRAL, elements of the quadrangular 
type are used exclusively because of their convenience. 
The computational domain is filled with a topologi- 
cally rectangular set of 

D-l 
q'=0 

of these building blocks, where D{< 4) is the dimen- 
sionality of the problem, and Hq,/ is the number of 
elements along each mesh dimension a' . Each ele- 
ment that borders the domain has one surface lying 
on the boundary that itself is an element of dimension 
D — 1. The total number of such boundary elements 
enclosing this rectangular mesh is a sum over the rect- 
angular faces 

D-l D-l 

^ = 2 E U^f^' (6) 

a'=Q P'^a' 

The element mesh can be distorted by stretching, 
compressing, bending, or even twisting it to conform 
to the geometry of the domain (as, for example, in 
a curvilinear coordinate system). In the engineer- 
ing FEM this coordinate transformation is called the 
"isoparametric" transformation, because coordinate 
values a;" and the variables are specified at the 
same nodal points. In general relativity this transfor- 



5 



mation is the generalized Lorentz transform 

fa _ 

with inverse 



where is the coordinate in mesh space, with range 
< C < 1 in each dimension o! . Basis vectors along 
the mesh coordinate direction a', and corresponding 
1-forms, are 

Appendix B discusses conditions that may need to be 
satisfied by this transformation. However, unless one 
wishes to use the mesh as an actual Lorentz frame of 
reference, or wants to follow the evolution of all wave 
phenomena, only the Jacobi condition is necessary for 
numerical stability 

C = det||£|| (9) 

2.3. The Choice of a Basic Coordinate Sys- 
tem 

In the past, when developing a finite difference nu- 
merical simulation code, for example, it has been cus- 
tomary (and considered necessary) to write the dif- 
ferential equations in the same coordinate system de- 
scribed by the computational mesh. That is, if the 
mesh is spherical-polar, then the equations are writ- 
ten in spherical-polar coordinates, and so on. How- 
ever, in the numerical method developed in this pa- 
per, this degeneracy is neither necessary nor desirable, 
as the mesh coordinate system is unknown until the 
problem is solved. 

To allow for an arbitrary, unknown mesh, the dif- 
ferential equations will be written in a "basic" or 
"real-space" (a;") system which does not change as 
the calculation proceeds. The derivatives still will be 
computed in the mesh ) system, but, in order to 
use them in the differential equations, will then be 
transformed to the basic system 

[VH"„ = W\a = C^'aW^o., (10) 

using the isoparametric/Lorentz transformation. 

The choice of coordinate system for the mesh is 
determined by how one lays out the elements in real 
space. That is, 

9a'l3' = J0."a' 9a0 (H) 



gives the metric coefficients in mesh space. Howc!vc;r, 
one still needs to choose a system in which to write 
the differential equations; but, since the computer will 
be doing all the curvilinear work for us, one can select 
a very simple basic system, keeping the coordinates 
as Cartesian (or as Minkowskian) as possible. For 
example, in axisymmetric problems, cylindrical coor- 
dinates will be used, not spherical-polar. For three- 
dimensional problems, Cartesian coordinates will be 
used, no matter how spherical the star or flattened 
the accretion disk. Any curvilinear properties of the 
metric orthogonal to the computational domain will 
be embodied in the volume element g, and curvi- 
linear behavior within the domain will be handled by 
the isoparametric transformation. 

2.4. Discrete Solution: The Nodal Mesh and 
Interpolation Scheme 

As with all continuum numerical methods, the so- 
lution is expressed as a finite set of discrete values. 
In the FDM these are values of the solution at speci- 
fied points (nodes) in space; in spectral methods these 
are coefficients of basis or interpolation functions. In 
the FEM, these discrete values are both nodal values 
and basis function coefficients. That is, the FEM has 
properties of both finite difference and spectral meth- 
ods. 

The finite element nodes are distributed within 

each (;lc;ment in such a way that the w;" can be in- 
terpolated across the element in each dimension with 
at least linear accuracy or better. For quadrangular 
elements, the simplest approach is to fill each element 
box with a (possibly hyper-) cubic mesh of {pa' + 1) 
nodes per dimension a', where pa' is the order of in- 
terpolation in that dimension, and nodes are shared 
by adjacent elements at all the interfaces (corners, 
edges, faces, and hyperfaces). For a problem of total 
number of dimensions D, the total number of nodes 
describing each element is, then, 

D-l 

I = U (Po' + 1) (12) 

a'=0 

That is, for four-dimensional elements. / = 16 for 
first-order (linear) interpolation, / = 81 for second- 
order (quadratic) interpolation, and / = 256 for third- 
order (cubic) interpolation — just within each ele- 
ment. The total number of nodes in the entire mesh 
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IS 



D-1 



2: = n (^"'po' +1) 



(13) 



(no sum on a'). The number of nodes on each ele- 
ment 's boundary is the total minus those in the inte- 
rior 



D-l 



K = I Y[ (Po' - 1) 



a'=0 



and for the entire mesh 



D-l 



jc = 1 - n(^«'p"'-i) 



(14) 



(15) 



No matter what the value of pa', in these simple 
cases the basis functions in mesh space for a node i 
within a given element e are products of Lagrange 
interpolation polynomials £iea' in each dimension a' 

_D-1 



£iea' (C ) 



a'=0 

D-l I I pa' _ pa' \ 

= nn|^ (16) 

where is the mesh coordinate value at node i in 
element e, AA'je is the contribution from that element 
to the basis or "shape" function for node ?, as mea- 
sured in the mesh (primed) system. The range of the 
nodal indices in the entire mesh is j = 1, ...,X, but 
the product in equation (|l^) only runs over the nodes 
within element. The total basis function for node i is, 
then, a sum over the element contributions 



(17) 



(which actually involves only those elements contain- 
ing that node). Note that each shape function attains 
unit value at its own node and zero at all other nodes 
in its associated elements (and in the mesh as well) 



(18) 



Often the body-centered nodes, and sometimes 
even face-centered nodes, are removed from the stan- 
dard Lagrangian elements elements to form the so- 
called "serendipitous" elements (Zienkiewicz 1977). 
In that case, if node t is removed, then the basis 



functions are given by the normal Lagrange shape 
function with the Lagrange shape function for that 
missing node subtracted off 



(19) 



with i ^ i. 

Most of the properties of Lagrangian elements can 
be illustrated in one dimension. Figure |^ shows a sim- 
ple 1-dimensional, 5-node mesh and its discretization 
in linear and quadratic elements. Note that interior 
shape functions have continuous derivatives at their 
respective nodes, while shape functions on element 
boundaries have discontinuous derivatives. (The lat- 
ter also involve more nodes as they are composed of 
shape function pieces from adjacent elements.) The 
FEM, therefore, can be considered to be a multi- 
domain spectral method with each of the thousands 
to millions of elements being a separate domain. 

Because the shape functions are continuous through- 
out il, the solution w and the coordinates x are truly 
continuous functions of position in the mesh: 



x"(0 



A^'?(0< 



(20) 
(21) 



like spectral methods but unlike the FDM where in- 
terpolation is only an ad hoc addition to the scheme. 
Also, because a unique inverse relation ^ = iix) ex- 
ists, the variables have an implicit function of position 
in real space 



(22) 



2.5. Formation of Derivatives and the Differ- 
ential Equations 

We now have a numerical procedure for computing 
the derivatives of the w'" with respect to x". First, 
the coordinate transformation matrix is formed 



/-a A/"/, 

a' — •'^ I, a' ■ 



(23) 



and then inverted to obtain Then the deriva- 

tives of the variables are computed in mesh space 



\a' = M'i^a'Wl 



(24) 



and, finally, transformed to real space by the chain 
rule 

W\a = C^' aM'l^a'Wl (25) 

Although not usually used in practice in the ac- 
tual computer code, it is sometimes useful for analytic 



7 



purposes to express the shape functions in real space 
coordinates and use them to interpolate the and 
compute their derivatives 



W'" ,a{x) = Mi,a{x)w^^ 



(26) 
(27) 



The real-space A/j also have the normalized property 
at their respective nodes 



By comparing (|2|) and (|27D with with (|2C|) and (|2^), 
one concludes that the real-space basis functions and 
their derivatives are 

Ui{x) = U'-Mx)) (28) 
UiAx) = C'^' cWi^c^'mx)) (29) 

With expressions for the w"" and w'" (either in 
mesh or real space) we now can calculate the residuals 
(equation [^) at any point in the domain f2 and 3?^ 
(equation ^) at any point on the boundary 957, not 
just at the nodes. 

2.6. Generation of the Nodal ("DifTerence") 

Equations: The Weighted Residual Method 

2.6.1. General Formulation 

The next step in the development of the astrophys- 
ical FEM is to construct a set of VI equations for 
the VI shape function coefRcients (wj ) that fully de- 
scribe w'"{x). This is accomplished by integrating 
the physical differential equations (|l]) and/or bound- 
ary conditions (^ over a function yVi(a;) which peaks 
near (but not necessarily at) node i and falls to zero 
far from that node. This produces VI discrete nodal 
equations 



"iqiiw)) = / Wi{x)^q{w])d^l = 



(30) 



n 



Relaxation schemes in the code then attempt to force 
the nonlinear 3gj to zero. In principle, each 3gj is 
a function of all of the w""-. However, in practice, 
because Wj is peaked near node S^j involves only 
nodes local to « — in fact, only nodes in those ele- 
ments containing node i. The 3^^, therefore, are more 
similar to difference equations than to spectral equa- 
tions and, when linearized, produce a banded rather 
than filled matrix. 



Because 3?^ can contain second derivatives, the in- 
tegral in equation (^0|) cannot be performed uniquely 
for every node using only the interpolation within a 
single element.^ The standard solution to this prob- 
lem, and the key step in the finite element process, is 
to integrate the weighted residual by parts to arrive 
at the so-called "weak" form 



an 



Wi{x)T^d{d^)p 



{Wi{x)^pTi; + Wi{x)Fq)dn 



(31) 



where, in the mesh system, the volume scalar is 



^"'^'T'^' W^"' An 



4! 



dC de dC dC 



(32) 



the surface 1-form normal to the domain boundary is 



= C^C-'p^2lEpLd^0' dC' df {S3) 



and e is the flat-space Levi-Civita permutation ten- 
sor. In the weak form, all terms involve only first- 
order derivatives of the variables with respect to 
the nodal coordinates. Second order derivatives are 
generated by the Wj^j term which, after integration, 
differences the flux on each side of each node in a 
manner similar to a finite volume scheme. 

Note that, in order to generate the nodal equations 
in the interior of the domain, weights yVj(a;) that van- 
ish on the boundary always will be used. Therefore, 
the first term in equation (31) — the boundary term 
— will always be zero. 

2.6.2. Second- order Equations: The Galerkin Method 

The weighted residual method can be derived in 
a number of ways. In early papers on finite element 
analysis, only linear problems were addressed and the 
nodal equations were generated using a variational 
approach that maximizes the norm of the solution 



•^The problem occurs at boundary nodes, where one needs infor- 
mation in adjacent elements. For example, in one dimension 
linear shape functions have zero second derivative, so they can- 
not represent the kernel at all. Quadratic shape functions do 
have a second derivative, but only one unique value within a 
given element, which is valid for the central node, but not for 
the boundary nodes. 
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dZienkiewicz 1977D . This led to the form (|30|) with 
the shape function itself as the weight 



W^{x) = Mix) 



(34) 



This choice for Wj is called the Galerkin method and 
is especially useful for second-order equations. For 
example, for the simple Poisson equation in one di- 
mension [w^xx ^ P = 0), with quadratic elements and 
uniform node spacing Ax, equation generates the 
following nodal equation at the central node i of each 
element 



4Ax Wi-i - 2wi + Wi+i pi-i + 8pi + pi+i 



10 



= 



which is similar to the finite difference form (wj-i — 
2wi+Wi+i) / Ax'^ — Pi = 0. (Somewhat more complex 
4th order difference equations are generated for nodes 
on element boundaries.) In general the Galerkin 
weighted residual method generates derivatives simi- 
lar to those expected in finite difference schemes (al- 
though in general geometry), but scalars and source 
terms are weighted averages of nodes surrounding i 
rather than evaluated exclusively at i. 

An important property of the shape functions hints 
at a more fundamental interpretation of the weighted 
residual method that is not discussed usually in the 
engineering literature. As the element volume SVl ap- 
proaches zero, the shape functions become good ap- 
proximations to the Dirac delta function 



lim Afi{x) 

sn^o 



ki 5il 5{x — xi) 



(35) 



(no sum on i) where fcj is a scaling constant of or- 
der unity, but generally different for each node i. 
Therefore, with the Galerkin method, the weights 
Wi in equation ( |30| ) are generalized approximations 
to 5{x — Xi) which, when integrated over a differen- 
tial equation, generate or "pick out" the correspond- 
ing difference equation near node i. As the number 
of finite elements approaches infinity, the 5gj more 
closely approximate the complete set of 3?^ defined at 
all points in O. Therefore, while originally derived 
for linear equations, the weighted residual merhod is 
valid for nonlinear problems as well. 

2.6.3. First-order Equations: Petrov- Galerkin Schemes 
and "Staggered Grids" 

The lack of a single, universally-applicable weight- 
ing function Wiix) is the main impediment against 



developing a truly general simulation code. One must 
always know the order of differential equation being 
integrated. For example, while the Galerkin scheme 
works well for second-order equations, it has the 
same pitfalls for first-order equations and fluid flow 
as centered-differencing schemes have in the FDM: 
leapfrogging, in which important terms in the equa- 
tions do not depend on variables at the node at 
which the integral ( pT| ) is evaluated {{dw/dx)i w 
{wi^i — ?«j_i)/2Ax), and two-point oscillations near 
shocks. 

Such problems can be addressed by using weight- 
ing functions other than the Mi. These are called 
Petrov-Galerkin schemes (Hughes 1987). For odd- 
order equations, functions that shift the peak of the 
weight away from node i reduce or eliminate many 
of these problems. This is the case for weights that 
are shape functions of twice the element interpolation 
order 

(36) 



PGl 



(x) = M^^^ix) 



and for those that are distorted by the shape function 
derivative 



wr\x) = M ± v-MiAx) 



(37) 



In the first case i + \ signifies a position in the mesh 
centered between nodes. The functions W/^*^^ peak 
in between nodes and have almost the same effect 
as using a staggered grid does in the FDM. Both 
W?'^^ and WP'^'^ generate non-leapfrogging differ- 
ences {{dw / dx)i w (wj+i — Wi) / Ax) and are still gen- 
erally second-order accurate (or higher) as the scalar 
source terms are evaluated at the same place as the 
derivatives. 

An upwinding scheme can be generated by setting 
v" (X u" and integrating only the advective terms 
with this weight. However, this scheme has l ow order 
accuracy and is rather diffusive ( Hughes 1987 ). Better 
methods for han dling shocks are the van Leer scheme 
( van Leer 1979 ), in which one enforces monotonicity 
in the gradient of the flux of a co nserved quantity, and 
highe r order Godon ov schemes ( Colella fc Woodward 



1984; Colella 1990 ) in which one applies the shock 



jump conditions within the element itself. Upwinding 
schemes in the FEM are a sub-field in themselves and 
are largely beyond the scope of this paper. 

2.7. Application of Boundary Conditions 



Not all of the nodal equations generated by (31) 
are useful. For example, for each second order differ- 
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ential equation (one where a given 1^ is a function 
of the gradient of at least one uf) exactly /C of the 
nodal equations (those integrated over shape func- 
tions peaking at boundary nodes on 90) are mean- 
ingless or incomplete. This is due to the absence of 
elements beyond the boundary needed to complete the 
integrals. For first-order equations (ones where T^^ is, 
at most, a function of w and x only) those on only 
one portion of the boundary must be discarded {e.g., 
on one side or surface). These ignored nodal equa- 
tions must be replaced by exactly the same number 
of boundary or initial conditions. These are gener- 
ated in a weighted residual manner similar to that in 
equation (|3^) except that the integral and weighting 
functions now are evaluated on the boundary 



dn 



Wi^{x)fr{w])d{dn) (38) 



where d{dU) is the magnitude of the surface element 
d{dn) = n^d{d9)p (39) 

and the nodes k lie on The number of equa- 

tions, therefore, will remain equal to the number of 
unknowns, as is necessary for a well-posed problem. 

2.8. Numerical Aspects of the Method 

2.8.1. Integration oj the Weighted Residuals 

In the FEM the integrals in equations (|l]) and ( |3^ ) 
often ar e performed numerically us ing Gaussian inte- 
gration ( Abramowitz fc Segun 1965 ), with all sampled 
points Xg interior to element boundaries. In prac- 
tice, the integrals are calculated piecewise, element 
by element, with each element's contribution to the 
various integrals summed accordingly. To accomplish 
this, one defines a local mesh coordinate system, refer- 
enced to the element center and parallel to the global 
^" system 



d 



— (; a" 



2(e - C ) 



-1 < sf < 1 



(40) 

(41) 
(42) 



where is the position in mesh space of element e's 
center, A^" is the element width in direction ^" , and 
double primes refer to the local element system. The 
Lagrangian shape functions within each element take 



on a very simple form in this local system. For linear 
interpolation in one dimension (with nodes lying at 
ST, = -1 and sa„ = +1) 



A/^'^(Se) = 0.5(1 + SieSe) 



(43) 



For quadratic interpolation (with element nodes lying 
at Sfe = —1, 0, +1) 



0.5 Sie Se (1 + Sie Sg) (sfe = ±1) 
(1 - 4) (Sfe = 0) 



(44) 

and so on for higher order interpolation. (Equations 
p3[ and Q are the functions depicted in Figure 0.) 
Higher dimensional element shape functions are prod- 
ucts of these in a manner similar to equations ( p^ ) for 
Lagrangian elements and (|l^) for serendipitous ele- 
ments. 

Equation ( ^l| ) then becomes, dropping the first 
term, as discussed earlier, and summing over elements 
and Gaussian integration points, 

"^H ~ - im^9 [Wi(a;ge),/5T,^ + WiiXge) Fq] Sfle 



where 



D-l 



(45) 



(46) 



a'=0 



is the element volume in the real-space coordinate sys- 
tem, e = 1, is the element number in the mesh, 
5 = 1, G is the number of the Gaussian integration 
point within element e, and ujg is the Gaussian weight 
at that point. Similarly, equation (p8[) becomes 



6^ 



- 5m ^9 Wfc i^ab) fr rib'^' 6{dnb)c.' (47) 
b g 

where 6 = 1, is the boundary element number, 

D-l 



rk 



(48) 



is the surface 1-form on the a' boimdaries, and is 
one of ^"'o'/V-5o'0', 5°"' v / \/-gvv , etc., depending 
on the boundary. 

Because the residual weights [i>Vi[xg) = yV"j(sg)], 
their derivatives, and the Gaussian weights ojg are the 
same for all elements, they can be precomputed and 
stored prior to beginning the relaxation of the solu- 
tion. The only quantities necessary to compute dur- 
ing the relaxation are the Tj^ and Fq at each Xg^ inte- 
rior integration point and the fr at each Xgb boundary 
integration point. 
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Historically, the number of integration points G 
used in each element is a function of the expected 
nonlinearity of the product WjSRq with respect to po- 
sition X. In engineering, this is usually of low order 
(linear or quadratic) , but in astrophysics this product 
can vary with exponential order or higher. Neverthe- 
less, in practice, even with highly nonlinear functions, 
the author has had quite satisfactory results using the 
same number of integration points as nodes in each 
dimension (i.e., G = I). Fewer than this ("underin- 
tegration" ) reduces the order of accuracy or even can 
produce a singular matrix. More than this ("overin- 
tegration" ) does little to improve accuracy (of order 
unity improvements only) at great computational ex- 
pense. 

2.8.2. Solution of the Simultaneous Nonlinear Equa- 
tions: The Multi- Dimensional Henyey Method 

For solving the VI nodal equations and bound- 
ary conditions, the author currently uses a standard 
multivariate Newton-Raphson technique, sometimes 
called the "Henyey" method in astr ophysics ( [Clay- 
ton 1968| ; [Aller fc McLaughlin 196^ ). The ^ai are 



linearized by expanding in a Taylor series about the 
solution , resulting in a matrix inversion problem 



(49) 



(no sum on the iteration number [n]) to solve for the 
corrections Swp"^^ which need to be appHed to the 
current guess to obtain a new [n + 1st] approximation 
to the solution 



["+1] 



Sw 



V [n] 



(50) 



In the engineering FEM d'^qi/dw'- is called the "tan- 
gent stiffness" matrix. It has a length VX on a side 
and bandwidth ~ VI^^-^/°\ At present GENRAL 
uses direct methods (Gaussian elimination with lower- 
upper decomposition) to solve equation (^), repeat- 
edly applying the corrections until the norm over all 
V and J 

A.„ = ||(5u;|["V«^i'"+''ll 
falls below a certain tolerance. 

The tangent stiffness matrix need not be extremely 
accurate. Indeed, when A^, << 1, the matrix need 
not be recomputed at all, with little impact on the 
rate of convergence and no impact on the accuracy of 



the solution. Furthermore, the elements of the stiff- 
ness matrix in equation ( ^9[ ) can be calculated using 
numerical differentiation, rather than writing out ex- 
plicitly the partial derivatives of each equation with 
respect to each variable. This eliminates the need 
to know the geometry of the mesh beforehand — a 
feature important for multidimensional astrophysical 
structures. Numerical differentiation of d^qi/dw'^ is 
not necessarily more expensive than algebraic differ- 
entiation, especially if the the baseline residual inte- 
grals 3^"' are calculated only once for each matrix, 
and those partial derivatives known to be identically 
zero are not computed. 

2.8.3. Logarithmic Variables 

It is quite common for an astrophysical state vari- 
able — e.g., the density p — to vary by many orders 
of magnitude over 17. Therefore, in order to maintain 
the same relative accuracy over the domain, it may 
be necessary to solve for a much more slowly varying 
function, e.g., p = log]^o(p)- addition, for variables 
that can be positive or negative (like velocity) one 
may need a more complex function and its inverse 



V = slogio(i;) 
sdex({;) 



S log^oil + Sv/Vscale) (51) 



S Vs 
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Sv 



1) 



(52) 



where S = sgn(i;) — sgn(z;) and Vg^aie is a fixed scaling 
value for v. These "scaled logarithmic variables" are 
linear for \v\ « v scale and logarithmic for \v\ » 
Vscaie and can be negative or positive. 

Unfortunately, there is some loss of convenience 
and intuition in re-writing the equations in terms of 
these new, modified variables, especially if there are 
many of this nature. Therefore, the following proce- 
dure has been devised in order that the differential 
equations still can be coded in their basic form (i.e., 
using p and v) while maintaining the accuracy of solv- 
ing for logio(p) and slogio(-y): 

1. Each variable is flagged as being linear, loga- 
rithmic, or scaled logarithmic, and then stored 



logio(w'') 
slogio(w^) 



(if linear) 

(if logarithmic) 

(if scaled logarithmic) 



11 



2. Partial derivatives of the integrals in equation 
( |49| ) are calculated with respect to w'". How- 
ever, when computing the functions , Fg, and 
fr, the stored variables and their gradients are 
re-exponentiated to their unmodified forms with 



w 



10™ 

sdex(tl;^) 



(if linear) 

(if logarithmic) 

(if scaled logarithmic) 



and 



With this scheme one can choose a variable to be loga- 
rithmic or not at runtime, or even switch its character 
during execution, without modifying the code. 

2.8.4. Pivoting 

Lower-upper decomposition techniques work well 
only when the matrix elements on the diagonal are 
decidedly non-zero. That is, one must identify which 
equation is "for" which variable w| . The term piv- 
oting refers to the exchange of rows and/or columns 
in the matrix to ensure that all elements on the di- 
agonal are indeed large — i.e., that the trace of the 
stiffness matrix is a maximum. 

Local, or partial, pivoting ensures that, at a given 
node, the correct physical equation is paired with the 
correct variable. For example, in MHD computations, 
when the electrical conductivity is infinite, the current 
jTis determined by Maxwell's equations, not Ohm's 
law. Or, in a hydrostatic star the momentum equa- 
tion determines the pressure structure, while velocity 
is determined by energy or particle conservation. In 
the standard partial pivoting algorithm one searches 
a matrix column for the largest element and switches 
the row of that element with the one presently on that 
column's diagonal ( Press et al. 1989| ). When apphed 
at a single node, this algorithm is very successful in 
automatically pairing equations and variables. That 
is, the solution will be able to evolve from a dynamic 
state to a hydrostatic one without re-casting the equa- 
tions or writing a new simulation code. 

Global pivoting ensures that an equation (in- 
tegrated near node i) is applied at the correct node 
J. This is a more difficult task than local pivoting 
and is not easily automated in our case. The correct 
identification is not always J = i, especially for first- 
order equations, for which the answer is determined 



by where the boundary conditions are applied. For- 
tunately, the global pivot of the matrix is a property 
that usually does not evolve with the simulation; it 
need be determined only once. Therefore, an elabo- 
rate pre-pivoting scheme has been devised for GEN- 
RAL, in which, based only on the equation order and 
location of boundary conditions, shape function co- 
efficients are identified with nodal equations, useless 
nodal equations are discarded, and boundary condi- 
tions are inserted. Local pivoting can still shift em- 
phasis for a given variable from one differential equa- 
tion to another, but the overall global identification 
of integrals (31) with nodes remains fixed. 



3. The Finite Element Method with Nodal 
Coordinates as Part of the Solution 

In most multi-dimensional astrophysics problems, 
the grid chosen initially will not be a good match for 
the final structure, due to a poor fit to the object 
boundary or poor resolution in areas of rapid gradi- 
ents. Therefore, the coordinates and/or quantity of 
nodes should be changed as the solution converges in 
order to get a better fit ("adaption"). In addition, 
one needs to allow for the mesh at one time step to 
be different from that at the previous time step (grid 
motion) . 

3.1. Adaptive Gridding 

Adaptive gridding, as used here, means modifying 
the mesh spacing in order to achieve greater accu- 
racy or stability without changing the total number of 
nodes or the topology of the mesh. One-dimensional 
stellar structure models use a form of adaptive grid- 
ding as they utilize a mass coordinate rather than 
radius, allowing the radius of each mass zone, includ- 
ing the outer stellar radius, to expand or contract de- 
pending on the current state of the star. Our general 
multidimensional adaptive gridding scheme takes the 
same form as equation (|l|) except that it is written in 
mesh space 

^"''^',/3' = (53) 

where y^" ^ is the adaptive gridding tensor. Cur- 
rently the author is using a diagonal form 



[/,(Ciea>-Vz/;"- ^ , )+- 



Ax" 



Ax"' 



(54) 



which ties the mesh spacing to the local gradient of 
the state variables. Note the sum over V variables; fy 
is a vector of I's and O's that, at run time, selects those 
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state variables to which the grid should be adapted, 
Ba' is a unit vector in spacetime along the local mesh 
direction 



ex.' 



(no sum on a'), <l^w^> is the average change in 
along the mesh direction ^" (determined from the 
current solution for w"), Ax" is a measure of the 
local linear mesh spacing along an element edge 



e-, ■ 

ISa'c'T/'Ar' 



(55) 



(again, no sum on a'), and Ci ~ 0.2 is a constant that 
regulates the strength of the gradient term. Equa- 
tions ( p3[ ) provide four constraints on the nodal coor- 
dinate values, allowing the mesh spacing along to 
decrease in regions of high gradients in the variables, 
but to be uniform otherwise. Note that the adap- 
tive gridding equations are not the coordinate condi- 
tions needed to complete the description of the metric 
(see Appendix A) . They simply move nodal positions 
around in an already-determined metric. C2 is a very 
small constant (« 10^^°), such that in the absence 
of adaptive gridding {i.e., all f" — 0) and with the 
proper boundary conditions, the mesh will assume an 
appropriate curvilinear character with uniform spac- 
ing Ax" in each . 

3.2. Moving Grids and the Advective Deriva- 
tive 

Adapting the grid to local conditions changes the 
spatial and time coordinates of the nodes. For an ob- 
server traveling along this produces what appears 
to be motion of the grid through space. Now, com- 
putational fluid dynamics on stationary, uniformly- 
spaced meshes is a fairly complex field in itself. On 
arbitrarily-spaced and moving grids the resulting mixed 
Lagrangian-Eulerian hydrodynamic equations might 
seem intractable. However, with the FEM the ex- 
act opposite is true. // the elements have an extent 
in time as well as in space, then the inverse isopara- 
metric transformation will automatically take into ac- 
count the complicated effects of differencing the fluid 
equations with respect to a moving grid. It is unneces- 
sary, and indeed incorrect, to attempt to include the 
grid velocity in the differential equations. 



As an example, consider the non-relativistic total 
advective derivative in flat spacetime 



D 
Dt 



— + v' — 

dt dx 



dx^ di° 



(56) 



where is the three-velocity of fluid flow. Now, let 
be a grid three-velocity such that e^;' remains 
parallel to e^i, although e^o' can make an angle 
(tan 6* 
tions 



At/|Ax|) with et- (Thes e are the condi- 



LeBlanc fc Wilson (1970, 1971) placed on their 



moving grid in their Lagrangian-Eulerian MHD calcu- 
lations.) The isoparametric transform and its inverse 
then are 



dxf_ 







AJ3" 

V' AC 



(if 13 ^ a') 

(if p^i;a' ^0) 
(otherwise) 

(if a' = 13) 



(57) 







Ak* 



(if a' ^ 0; P = i = i')^^^'^ 
(otherwise) 



(no sum on (3 or i) 
(p6|) one obtains 



Substituting equation (p8|) into 



D 
Dt 



AC d 
Ax* dC' 



(59) 



The first term in equation ( |59| ) is the apparent time 
derivative (along ) at a given node in the grid 
frame while the second is simply the moving advec- 
tive derivative {v — Vg)-'V. Thus the isoparametric 
transform reproduces the Lagrangian-Eulerian equa- 
tions of LeBlanc & Wilson under the same conditions. 

3.3. Adaptive Mesh Refinement 

The purpose of adaptive mesh refinement (AMR) 
is similar to that of adaptive gridding: to increase 
resolution in a local region of the mesh. However, in 
the case of AMR the numbers of elements and nodes 
usually increase as the calculation proceeds. When 
the element mesh extends into the time dimension, 
AMR makes it possible to treat eflticiently the prop- 
agation of particularly important wave phenomena. 
Rather than subdivide the entire domain finely, one 
does so only for the region of spacetime near the null 
geodesies along which the phenomenon propagates. 
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The Courant condition (Appendix B) prescribes one 
form of AMR, as it places an upper limit on the tem- 
poral mesh spacing and, therefore, a lower limit on 
the number of elements necessary to solve a problem. 
If the Courant condition is not satisfied in a local re- 
gion, the number of elements in ^'^ will have to be 
increased. 

In the FEM AMR can be achieved by subdivid- 
ing some of the elements, usually by a power of two 
in a given dimension, and then rebuilding the inter- 
polation grid of nodes within the newly-gridded sub- 
domain. There are two matters of concern in this 
case: 1) the subdivision must be done by elements, 
not nodes, regardless of the final order of the inter- 
polation grid used within each element, and 2) larger 
elements bounding the more finely subdivided region 
must have their internal interpolation scheme, and 
hence their nodal shape functions, modified by the 
addition of new boundary nodes. That is, these mod- 
ified elements must be serendipitous elements. This 
ensures that the functions w'"{^) and are con- 

tinuous across element boundaries and that there are 
no spatial gaps or "hanging nodes" . This second re- 
quirement can significantly increase the complexity of 
the mesh and the number of element types for which 
quantities need to be pre-computed and stored. 

Recently, however, the astrophysical community 
has been embracing an alternative, multi-^et;e/ ap- 



proach to AMR ( [Truelove et al. 1998| ; [Norman 19981 ). 
A region of space is refined not by subdividing the 
grid cells themselves, but by applying separate, and 
successively more refined, grids at the same location 
and with some nodes in common between each level. 
This hierarchical approach eliminates the hanging 
node problem without resorting to defining many new 
serendipitous element types: each grid level is sub- 
divided with standard linear or quadratic elements. 
This technique also fits naturally into a multigrid iter- 
ative scheme for solving the coupled nodal equations. 

4. Tests and Examples 



Below are presented some tests of the code GEN- 
RAL on problems with known solutions. At the 
present time the code uses the multi-dimensional 
Henyey technique to solve the difference equations it 
generates. As discussed earlier, this approach has se- 
vere computer time and memory limitations. There- 
fore, all the examples involve a much smaller num- 
ber of elements than the millions that one would use 



in a typical astrophysical simulation. (Generally, the 
full code running on a desktop workstation is lim- 
ited to 4096 elements or 6561 nodes total [64^, 16^, 
or 8** linear elements or 32^, 8'^, or 4^ quadratic ele- 
ments]). Nevertheless, the tests serve to demonstrate 
the unique features of the FEM, including its ability 
to solve nonlinear astrophysics- like problems in multi- 
dimensional, arbitrary curvilinear coordinate systems 
and to achieve high accuracy in the solution by em- 
ploying higher order interpolation, adaptive gridding, 
and logarithmic variables. 

4.1. Simple Tests 

J^.l.l. Tests of the Interpolation Scheme 

The first tests of the method involve no differential 



equations. Instead, the integrands in equations (45) 
and ( ^7| ) are set to unity so that only domain volume 
and surface area are computed. The results are then 
compared with known solutions for various shapes 
(rectangular hyper-solids, hyper-spheres, and ellip- 
soids), different element interpolation schemes (lin- 
ear [p = 1] and quadratic [p = 2]), different num- 
bers of elements H in each dimension (1,2,4,8,16), and 
different numbers of integration points G (p^ and 
(p -I- 1)^). This procedure is more than simply a 
verification that the numerical integrations work. It 
also exercises the coordinate transformations and the 
ability to form proper volume and surface elements in 
fairly arbitrary curvilinear situations. 

Table 1 summarizes the results, normalizing to 
grids with nine nodes in each dimension.^ In this case 
(but, unfortunately, not for the differential equation 
tests below) both choices for G produce the same re- 
sults: for linear elements the error Eji = (ilf< — Vt)/Vl 
varies as Ax^ cx (Hp — 1)~^, while quadratic ele- 
ments produce fourth-order errors (Eq cx Ax'* cx 
(Hp — 1)~*). Results for surface area integration are 
similar, although the actual fractional error is about 
a factor of two smaller. Results for other ellipsoidal 
objects are also similar. 



''Given that the method is element-based, one could compare 
equal numbers of linear and quadratic elements in each dimen- 
sion. In that case, the errors in the right-hand column of Table 
1 would be a factor of 16 smaller. However, performing the 
comparison in the latter manner would make it difficult to tell 
which part of the improvement for quadratic elements is due to 
the interpolation scheme alone and which is due to the increase 
in the number of nodes. 
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TABLE 1 

Error Results for Volume Integration of D-Spheres'' 
(Meshes with 9 Nodes in Each Dimension) 



D 



Fractional Error for Element Type 
Linear*^ Quadratic^ 



2 
3 
4 



-6.4 X 10^3 
-2.0 X 10-2 
-3.7 X 10-2 



-4.9 X 10-'^ 
-9.0 X lO-'^ 
-2.2 X 10-4 



"Surface area results have a similar scaling, but with a factor of 
two better accuracy 

''Accuracy for linear elements varies with Ax^ and quadratic 
elements as Ai* 



TABLE 2 

Error Results for Cartesian Grid Tests 
(Meshes with 9 Nodes in Each Dimension) 

Normalized L2 Error for Element Type 
n D Linear*^ Quadratic^ 



1 


1 






b 




b 




1 


2 


2.5 


X 


10- 


-3 


1.2 X 10" 


-3 


1 


3 


1.4 


X 


10" 


-3 


4.0 X 10- 


-4 


1 


4 


9.9 


X 


10- 


-4 


1.8 X 10- 


-4 


2 


1 


5.8 


X 


10- 


-3 


b 




2 


2 


3.3 


X 


10- 


-3 


b 




2 


3 


2.3 


X 


10- 


-3 


b 




2 


4 


1.8 


X 


10- 


-3 


b 




3 


1 


1.2 


X 


10- 


-2 


1.2 X 10" 


-3 


3 


2 


6.1 


X 


10- 


-3 


6.2 X 10- 


-4 


3 


3 


4.0 


X 


10- 


-3 


4.0 X 10- 


-4 


3 


4 


3.0 


X 


10- 


-3 


2.9 X 10- 


-4 


4 


1 


2.1 


X 


10- 


-2 


3.1 X 10- 


-3 


4 


2 


1.0 


X 


10- 


-2 


1.6 X 10" 


-3 


4 


3 


6.4 


X 


10- 


-3 


1.0 X 10" 


-3 


4 


4 


4.6 


X 


10- 


-3 


7.2 X 10- 


-4 



"Accuracy for linear elements varies with Ax^ and quadratic 
elements as Ax^ 

''Exact solution obtained for one-dimensional first-order prob- 
lems and for all second-order problems with quadratic elements 



4.1.2. Fixed Cartesian Grid Tests: Poisson's Equa- 
tion in One to Four Dimensions 

The second set of tests involves solving a differen- 
tial equation in up to four dimensions, but on a reg- 
ular Cartesian (not Minkowskian) grid of dimension 
D. The equation used is Poisson's equation 



with a known solution 

Wo = r" 



(60) 



(61) 



where r = x"^ + y'^ + + t'^ is the radial distance 
from the origin and < < 1. The components of 
the generalized equation dll) are then 



T'^ = w^f} (62) 
F = V^wo ^ n (n + L> - 2) r("-2) (^53) 

with Dirichlet conditions on the entire boundary 

f = w- [r{dn)]" (64) 

where r{dfl) is the expression for r evaluated on the 
boundary. This test exercises the code's ability to 
solve equations ( p5| ) and (p7[), but does not test the 
coordinate transformations. 

Results of the fixed grid tests are given in Table 
2, which shows how the accuracy of the solution in 
different dimensions, determined by the normalized 
"L2" error norm 



zL2 



Wo) 



1/2 



1/2 



varies with number of nodes and elements used. Two 
differences from the volume and surface area integra- 
tion tests are worth noting. Firstly, the solution er- 
rors for quadratic elements are third-order accurate 
(E^2 QT _ 1)^3-)^ j^Q^ fourth-order. Sec- 

ondly, underintegration (G = p^) does not work. G 
always must be equal to or greater than (p+ 1)^ (i.e., 
> 2^ for linear elements and > 3^ for quadratic el- 
ements) in order to generate the proper second-order 
differences in the Laplacian. Underintegration, at 
best, reduces the accuracy of the solution by one or- 
der. At worst, it destroys nearest neighbor differenc- 
ing, producing leap-frogged differences, and can lead 
to no solution at all. "Iso-integration" (G = /) is 
probably sufficient for any equation of the form (^, 
but it should be checked in each circumstance to be 
certain. 
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4.2. Adaptive and Curvilinear Grid Tests 

The third set of tests exercises nearly ah the fea- 
tures of the code in order to obtain a solution to 
a rather pathological Poisson problem — a Fermi- 
Dirac-like function 



1 



Wq 



(65) 



with a cold temperature of 1/f = 0.02. Such sudden 
exponential drops in the solution are common at stel- 
lar core-halo boundaries or stellar surfaces and are 
difficult to resolve accurately without a large number 
of nodes or a variable change (to optical depth, for ex- 
ample) . For this demonstration the conservative form 
for the stress and force terms (in up to four Cartesian 
dimensions) has been chosen 
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(no sum on /3) where a/3 
and 



(l,a, 6, c) are constants 



(For example, in one dimension, a = = c = 0; in 
two dimensions, b ^ c — 0; and so on.) And, as 
one is still interested at this stage in testing the FEM 
machinery and not the astrophysical viability of the 
code, once again simple Dirichlet boundary conditions 
are employed, so equation (^ becomes 
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exp[/(r(917) - 0.5)] 
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(68) 



rather than, for example, a multipole expansion of the 
interior solution. The above conservative form (K 
was chosen in favor of other forms (such as 
I) because its solutions are particularly accurate 
for a small number of nodes and suitable for demon- 
strating adaptive gridding and the use of logarithmic 
variables. 

Figure || shows the solution of this Poisson prob- 
lem in one dimension as one applies successively more 
features of the code. The top panels of Figures || show 
standard fixed, equally-space grids of 8 and 16 linear 
elements respectively. Some improvement in accuracy 
can be obtained by doubling the resolution, but this 
incurs additional storage and computational expense. 



Turning on the adaptive gridding equations (q3|) , how- 
ever (middle panels), significantly improves accuracy 
for the same number of elements. This also demon- 
strates one aspect of the ispoarametric transforma- 
tion: variable node spacing. A closer examination of 
this more accurate solution, however (Figure ^, bot- 
tom left), shows very large relative errors in the log 
when w << 1. Nevertheless, these can be overcome 
easily, without re- writing and re-coding the equations, 
as the bottom right panel of Figure || shows. When 
w is identified as a logarithmic variable, rather than 
linear, the solution remains accurate over eleven or- 
ders of magnitude. (Note the different node spacing 
in the bottom panels, with the grid adapting to w on 
the left and w [— log^o w] on the right.) 

Figure ^ shows a similar development for a two- 
dimensional Poisson problem where the Fermi-Dirac 
surface has an elliptical ratio of 2:1 (a = 2). The 
first three panels demonstrate the errors possible in 
locating the surface if the proper coordinate geometry 
is not used. Additional improvement in accuracy can 
be obtained by using adaptive gridding (middle right 
panel). However, this solution for r > 0.5 suffers the 
same oscillatory errors seen in the one-dimensional 
case (bottom left), which again can be eliminated by 
identifying w as logarithmic (bottom right). 

It is important to note that, in all of the solutions 
displayed in Figures ^ and ^, no explicit curvilinear 
coordinate system is used. The coordinates of the grid 
points (whether fixed or part of an adaptive gridding 
solution) are stored only as xj and t/j, not as rj and 
9i, for example, and yet are still fully arbitrary (sub- 
ject to the Jacobi condition). The Poisson equation is 
written only in terms of a; and y as well. Of course, the 
derivatives are still calculated using the coordinate 
grids shown, but then they are immediately trans- 
formed to {x, y)-space using the transformations (^ 
and (|^). Thus, with these new techniques the grid 
can be moved around to obtain a more accurate solu- 
tion while the physical equations remain coded in the 
same very simple form. 

The Fermi-Dirac Poisson tests were used to deter- 
mine a good value for the adaptive gridding constant 
in equation (^). Several dozen models like those in 
Figures |2| and 3 were computed for different values 
this parameter. It was found that the accuracy im- 
proved by factors of 3 — 10 as Ci was increased from 
to 0.2, but beyond this point the accuracy did not im- 
prove much. In fact, for values much greater than 0.2, 
the models became unstable, often not converging. 



16 



Therefore, Ci — 0.2 was chosen as a semi-universal 
value in the adaptive gridding equation. It has proven 
to be useful both in the Fermi-Dirac tests in Figures 
H and 1^ and in the stellar structure models below. 

One important point about adaptive gridding should 
be mentioned. As currently implemented, the tech- 
nique is rather volatile and unstable. Unless great 
care is taken, iterations with adaptive grids often di- 
verge, violating the Jacobi condition in the process. 
In the case where a single solution to a steady state 
problem is sought, sometimes less CPU time cost will 
be incurred by subdividing the mesh more finely or 
using quadratic elements, rather than using adaptive 
gridding techniques. On the other hand, when many 
thousands of successive models are to be computed, 
as is the case for evolutionary problems, each newly- 
converged model will be a good initial approximation 
to the next evolutionary state, yielding convergence 
for each time step in only a few iterations. In this 
case, the amount of time spent converging the first 
adaptively-gridded model will be a small cost com- 
pared to the CPU time adaptive gridding saves over 
the course of the evolution by using a smaller number 
of elements and nodes to obtain the same high level 
of accuracy. 

4.3. Stellar Structure Tests: Polytropic Stars 

The fourth series of tests adds the ability to solve 
a coupled set of both first- and second-order partial 
differential equations. It also demonstrates the use 
of ^/—g to solve a problem in which the basic coor- 
dinate system (not just the grid) is curvilinear, due 
to the assumption of symmetry conditions involving a 
coordinate direction orthogonal to the computational 
domain. 

4-3.1. Spherical Polytropes in One Dimension 

In one dimension the equations for polytropic stel- 
lar structure are hydrostatic equilibrium (Euler's equa- 
tion with zero velocity), Poisson's equation for grav- 
ity, and the polytropic equation of state 



dv , , , dw 
-f + p{n + l) — 
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l+l/n ^7^) 



Pressure p and density p are unity at the stellar center 
and zero at its surface, and r is the spherical radius 



coordinate. The polytropic index n is a measure of the 
hardness of the equation of state; the factor (n -I- 1) 
in the hydrostatic equilibrium (HSE) equation is a 
normalization constant for the gravitational potential 
w. 

In semi-analytic treatments, the hydrostatic equi- 
librium equation is multiplied by r^/p, differentiated 
with respect to the radius r, and combined with Pois- 
son's equation to give a single second-order equation. 
Following this procedure here, however, tests only our 
ability to solve that simple equation and not much 
else. A slightly stronger test of the method would be 
to leave the system as a set of coupled equations and 
identify 
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P,r + p{n + l) 




: W^r (72) 

^ P (73) 
= (74) 

with boundary conditions at r = of p = 1 and 
p,r ~ 0. Note the need for a non-unit volume element 
of in order to form the proper divergence. A basic 
curvilinear coordinate system must be used because 
of the spherical symmetry assumed in directions or- 
thogonal to Br. 

Unfortunately, the above approach is still unsat- 
isfactory, because it cleverly casts HSE as a second- 
order equation, avoiding the first-order equation prob- 
lems discussed earlier. Also, as shown below, it does 
not lend itself to generalization to two or more di- 
mensions. To address these issues, the following al- 
ternate identification of the stress and force terms for 
the pressure equation has been chosen 



rpr 

p ~ 







Fp ~ p,r + p(n+l) W^r 



(75) 
(76) 



which casts HSE as a first-order equation (with only 
the boundary condition p = 1 at r = 0). Boundary 
conditions on the potential w are w^r = at the stellar 
center and w = —rsurface/r at the stellar surface. 
In addition, the condition p = ps = 10"** is applied 
at the surface on the adaptive gridding equation to 
determine the stellar radius. 

Figure ^ shows an n = 1 polytrope in one dimen- 
sion, which has the analytic solution 



w = —p — 1 



P = P 
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The left panel uses only the Galerkin method, while 
the right two panels show the solution using the 
Petrov-Galerkin scheme in equation (|3^) to integrate 
the HSE equation. The need to treat first-order equa- 
tions differently from second-order ones is clearly ev- 
ident. The leapfrogging first-order differences pro- 
duced by the Galerkin method not only display point- 
to-point oscillations, they also miss an implicit bound- 
ary condition {p^r = at r = 0, implied by w^r = 
and equation ^6|). The Petrov-Galerkin scheme im- 
proves the accuracy considerably, and with adaptive 
gridding remains roughly second order accurate for 
linear elements and third order for quadratic. 

At the present time no robust, automated method 
for determining the order of the differential equation 
has been developed. The code must be told explicitly 
not only the order of each equation but also the loca- 
tion(s) of the boundary conditions. While this is the 
greatest obstacle to producing a truly general contin- 
uum simulation code to solve all types of equations, 
it is a relatively modest amount of effort compared 
with writing a new code for each problem. 

4-. 3. 2. Rotating Polytropes in Two Dimensions 

Treatment of a uniformly-rotating polytropic star 
in two dimensions poses several problems. Firstly, it 
appears to be an overdetermined system 
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(where -Ejot = lu'^ /AnGp is the normalized rotational 
energy per unit mass, uj is the uniform angular veloc- 
ity of stellar rotation, and R and Z are the cylindrical 
radius and axis coordinates) with four equations but 
only three unknowns. In normal astrophysical situa- 
tions this dilemma will not arise, as the fluid equa- 
tions plus the conservation of energy are a well-posed 
problem. However, it occurs here because two un- 
knowns {vr and vz) and only one equation (conser- 
vation of mass) have been removed from the full set. 
The trick is to convert the two redundant equations 
for HSE (|^ into only one. 

One possible solution is to form a single second- 
order equation by taking the divergence of the HSE 
equation, similar to the standard semi-analytic ap- 
proach. However, while it works fine in one dimen- 
sion, this approach is unstable to two-dimensional 



perturbations when both the Dirichlet and Neumann 
boundary conditions are applied on the same (r = 0) 
surface, leaving the stellar surface free. Settings — pg 
at the surface does not help either; in this case the 
mesh must become adaptive, and this constraint must 
be used as a boundary condition on the adaptive grid- 
ding equations, not on HSE. Moving the Neumann 
condition to the stellar surface is a better approach, 
but difficult to apply for more complex problems {e.g., 
rotating polytropes). 

An approach that does work is to project the HSE 
equation along a direction in which p and w have 
significant gradients. The projection direction need 
not be along Bp = — Vp/|Vp|; appears sufRcient, 
even when the polytrope is rotating rapidly. But it 
must not be orthogonal to Bp, along which the gradi- 
ents are zero. The components of the general stress- 
force equation for the rotating polytrope are, there- 
fore, similar to (|7|)-(|7^) and (f|)-(||) 
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with a sum on i over R and Z, boundary conditions 
p = 1 and e^w^i = at the stellar center and at the 
stellar surface 
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where << 1 is a small fraction of the central pres- 
sure and Ws is the specified surface potential. 

Two different methods were tested for calculating 
Wg. The first was an exterior multipole expansion 
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1=0 



where the are Legendre polynomials and Mi are 
the moments of the mass distribution in the star 



Ml = 



p{R,Z) r'^ RdRdZ 



which, because of additional equatorial plane symme- 
try, are non-zero only for even Generally, orders 
up to L = 12 were used. The second method used a 
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Green's function integral over the domain 



P 



n \x ~ Xs\ 

p{R,Z) RdRdZ 



X L, In 



i?2 + (Z-Z,)2]l/2 



iR-Rsr + {z-z,y 



(89) 



This expression is vahd for all continuous, Newtonian 
self-gravitating, axisymmetric systems.|^ The single- 
parameter complete elliptic integral 



ds 



(90) 



(1 — tanha cos27rs)i/2 

represents the summed relative contributions to the 
surface potential at {Rs,Zs) from different angular 
elements of a ring of matter at {R,Z). Ig diverges 
with a (a measure of how close the ring is to Xg), 
but it can be evaluated numerically easily and tabu- 
lated to a part in 10^° accuracy for the useful range 
< a < 10 {i.e., for rings approaching within only a 
fraction 2 x e~^° = 9.1 x 10^^ of |a;s|) over which the 
integral lies in the range 1.0 < le < 5.1. In all mod- 
els presented here, even with the largest meshes (33^) 
and adaptive gridding {Rs/SRs ^ 200), a remains 
well below 8.0 and le below 4.2. 

As implemented in the author's code, the speed 
of the integral outer boundary condition technique 
was significantly slower than the exterior multipole 
expansion, increasing the time to form the stiffness 
matrix (although not affecting the time to invert it) 
by factors of several. While the author made no 
attempt to optimize the implementation, even after 
such efforts it nevertheless should remain somewhat 
expensive, as it requires (for each non-zero stiffness 
matrix element and each right-hand-side vector ele- 
ment) the complete integration of the potential (equa- 
tion |89|) at 4-6 surface points, each integration being 
the equivalent of of a multipole moment computa- 
tion. However, while the multipole expansion began 
to break down for modest rotation speeds, the inte- 
gral technique converged with no problem for all rota- 
tion speeds up to breakup, making the extra compu- 
tational effort worthwhile and necessary. In addition. 



^The integral in equation ( pE| ) is over both northern and south- 
ern hemispheres of the star. If, as is often the case, reflec- 
tion symmetry is assumed at the equator and Q refers only to 
the northern hemisphere, then the integral must take into ac- 
count the southern contribution as well by summing one term 
as above plus one term with (Z — Zs) replaced by (Z + Zs). 



as the integrals are done only for surface points, this 
hybrid technique (differential equation with integral 
boundary conditions) still will be much cheaper than 
computing the global potential by performing such an 
integral for every point in the domain. 

Figure [s] shows two dimensional n = 1 non-rotating 
polytropes for the two element classes each with 9^ 
nodes — the analog of the middle and right panels 
of Figure ^ — using the multipole boundary condi- 
tion. Note especially the variable grid spacing near 
the stellar surface and the difference in smoothness 
between linear and quadratic interpolation. Figure |^ 
shows the same models for the 33^ node cases. Note 
especially the departure from sphericity in the pres- 
sure contours in the models employing linear elements 
that is absent in the quadratic element cases; the 9^ 
quadratic model is more accurate than the 33^ linear 
model. However, errors in the two cases scale only as 



^quadratic ^ (^p) 



which is one full order less accurate than expected. 
This may be due to the reflective boundary condi- 
tions along the axis and equator, which are only first 
and second order accurate, respectively, in the linear 
and quadratic cases. The nature of this boundary 
condition cannot be improved at this time, but will 
be once iterative/niultigrid methods for solving the 
coupled equations are implemented. This will allow 
boundary "ghost" elements to be handled easily, al- 
lowing application of boundary conditions as accurate 
as the mesh interior. 

The Maclaurin spheroid sequence (a series of n = 
0, uniform density polytropes) provides a full two- 
dimensional test of the method. (Note that the 
method takes no advantage of the uniform rotation, 
uniform density, or polytropicity, so the ability of the 
code to solve the Maclaurin problem is a good indica- 
tion of how it will do on general multidimensional stel- 
lar structure and other more complex problems.) The 
appearance of Maclaurin spheroids is similar to that 
of the n = 1 polytropes in Figure || and 6, but the log- 
arithm of the pressure varies little with radius except 
near the surface, where the radial grid spacing de- 
creases dramatically due to the sudden pressure drop. 
For this reason it is sufficient to use a larger pressure 
boundary value (ps ~ 10^^ rather than 10^^) in or- 
der to determine the locus of the Maclaurin spheroid 
surface. Figure |^ shows the analogy of Figure || 
for a Maclaurin spheroid with uo"^ /2-KGp = 0.224 — 
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very near the theoretical hmit of 0.2246656 (Tassoul 
(1978)). The integral boundary condition ( ^9| ) is used 
to compute the surface potential. 

The complete Maclaurin series from [x)^/27rG/3 = 
to 0.224 was computed in four different ways, using 9^ 
and 17^ nodes for the two classes of elements (linear 
and quadratic). Values of r = Erotl\Ew\, the ro- 
tational flattening ratio of the semi-minor and semi- 
major axes f ~ 1 ~ az/aR, and total angular mo- 

mentum J = E^^^ p(R, Z)R^ dil have been com- 
puted from these models and are compared in Figure 
~] with analytic curves from [Tassoul (1978) and Chan- 



dra sekhar (1969). The fractional errors for the four 
series are shown in Figure |^. The models are quite 
accurate for such a small number of elements, with 
errors in the 10"'^ to 10^^ range in the 17^ quadratic 
case. They show the expected result that third order 
interpolation is significantly more accurate than sec- 
ond order, but the decrease in the errors with increas- 
ing numbers of elements is not as steep as expected. 
Most of this behavior is probably due to the less accu- 
rate reflective boundary conditions mentioned earlier. 

5. Discussion 
5.1. Summary 

This paper has developed a general method for 
solving multidimensional structural, and dynamical, 
problems of astrophysics. Virtually all situations in- 
volving continuous media are potentially addressable 
— in normal flat Cartesian space or in curved space- 
time. Problems in this area include, but are not lim- 
ited to, the full structure and secular evolution of vis- 
cous, rotating (and even magnetized) stars and ac- 
cretion disks in two and three dimensions, interact- 
ing binaries, asymmetric stellar envelopes and winds, 
non-radial pulsating stars, nonlinear development of 
secular and thermal accretion disk instabilities, and 
stationary or evolving spacetimes. 

While this method is most useful for structures 
that evolve on timescales long compared to a dynam- 
ical time, there is no formal restriction on how short 
the evolution time must be. Therefore, the approach 
to dynamical instabilities from a stable conflguration, 
and even initial dynamical development, also can be 
studied, although the author still recommends the use 
of an explicit code for full dynamical evolution. 

The equations of continuum astrophysics have been 
condensed into a general compact covariant form, and 



that form encoded into the author's FEM program 
GENRAL. A user can solve a particular astrophysics 
problem by supplying a single subroutine that takes 
as input one given coordinate position, plus the value 
of the variables and their gradients there, and returns 
as output the differential etinaiions (i.e., the general- 
ized "stress tensor" , "body force vector" , and possible 
boundary conditions appropriate for that problem). 
The program then generates the nodal or "difference" 
equations on a user-specified general curvilinear grid 
using the finite element weighted residual integrals, 
and solves the large set of coupled equations to pro- 
duce the solution to the equations. While described 
by discrete nodal values, as in finite difference meth- 
ods, the finite element solution is continuous, as in 
spectral methods, due to the finite element interpo- 
lation functions. Either second order (linear interpo- 
lation) or third order (quadratic interpolation) accu- 
rate solutions are possible in the code. In addition, 
the positions of the nodes themselves can be part of 
the solution (a "rubber" mesh), allowing grids to be 
fit to unknown boundary shapes and regions of high 
gradients to be more finely resolved with the same 
number of mesh points. 

While the method is cast in a full covariant form, 
it is anticipated that the initial applications will be 
mainly in the area of non-relativistic stars or accre- 
tion disks in static gravitational fields. The covariant 
form, however, is important even for non-relativistic 
problems. When the mesh extends into the time do- 
main, even only for one or two elements, the coor- 
dinate transformations that are a natural component 
of the finite element method automatically generate 
any arbitrary Lagrangian-Eulerian (ALE) advective 
derivatives needed to take possible grid motion into 
account. 

The method has been demonstrated on astrophysi- 
cally interesting problems (spherical or rotating poly- 
tropic stars) in one and two dimensions, with full 
adaptive gridding, and on simpler problems in three 
and four dimensions. 

5.2. A Note on the Solution of Elliptic Po- 
tential Problems 

A great deal has been written on the numerical so- 
lution of astrophysical potential problems like equa- 
tion d?^). The technique used here has elements of 
past approaches plus some new features and is well- 



suited to the FEM. Like many past authors (Clement 



1978; Bonazzola, Gourgoulhon, & Marck 1998), the 
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approach here casts the problem as a differential equa- 
tion with boundary conditions specified on the exte- 
rior of the domain. However, rather than being a 
simple 1 jr or low-order multipole potential at a large 
radius in the vacuum region, the author's preferred 
boundary condition is an integral solution of the dif- 
ferential equation, specified at the stellar surface. 
This integral is physically equivalent (on that surface) 
to the "full integral" technique that computes the po- 
tential throughout the computational domain using a 



ferential equation itself ( 


rassoul (1978) 




Eriguchi & 


Mii 


Her 1985; 


Komatsu, Eriguchi, & Hachisu 1989 


)■ 



However, the author's method of evaluating this inte- 
gral is somewhat different as it requires no expansion 
in terms of Legendre polynomials, relying instead on 
a single, slowly- varying, numerically-tabulated func- 
tion le (a) to handle the axisymmetry of the potential 
field. When calculated in this manner, using the inte- 
gration techniques already available in the FEM, the 
numerical integral is a solution of the discrete finite el- 
ement equations themselves to within the truncation 
error. The boundary condition and interior differen- 
tial equation, therefore, match well, leading to good 
convergence of the models. 

This technique will work in any situation where the 
full integral technique can be used: extension of the 
Newtonian case into three dimensions will be trivial, 
and it will be straightforward for the general rela- 
tivistic case as well. In three dimensions there is no 
axisymmetry, so Ie{a) will not be needed, and |a; — a^gj 
will take the simple Pythagorean form. For axisym- 
metric relativistic stars, the four metric potentials are 
given by three Green's function integrals plus a first 



order equation (Komatsu, Eriguchi, & Hachisu 1989) 



Therefore, the three elliptic equations can be solved 
in the same manner as Poisson's equation is solved 
here (although probably using different tabulated 1^ 
functions), preserving the differential equations in the 
stellar interior but using the integral solution on the 
surface. The fourth equation would be handled with 
the Petrov-Galerkin scheme demonstrated in section 
4. The advantage of this approach compared to the 
full integral technique is speed (the number of vol- 
ume integrals is proportional to the domain surface 
area, not the volume) . The advantage over the multi- 
domain technique is convenience (one does not have 
to deal with the vacuum region and the matching of 
stellar and vacuum solutions). 



5.3. Unresolved Issues 

While the code and method are mature enough 
to begin solving two-dimensional structural problems 
routinely, there are several unresolved issues, men- 
tioned in the text, that must be addressed more com- 
pletely before the full potential of the astrophysical 
finite element method can be realized. 

First and foremost are the execution speed and 
memory issues. While the reader may consider the 
generation of the transformations and finite element 
integrals rather time-consuming, by far the greatest 
use of computer resources is the technique currently 
used to solve the coupled equations — the Henyey 
technique. For large meshes in three or more dimen- 
sions, it becomes prohibitively expensive, requiring 
thousands to many millions of years of CPU time 
(cx [Hp]'^'^^^) and equally absurd amounts of mem- 
ory (oc [Hp]^'^^^) to invert once. However, multigrid 



methods (Brandt 1977) need only about twice the 
grid size in storage (oc [Hp]^) and only require a few 
sweeps of the mesh to converge (c>c [Hp]^ log[Hp]^). 
The author and P. Godon have been experimenting 
with modern parallel multigrid algorithms in finite 
difference codes with considerable success. The CPU- 
time and memory scalings, and linear speedup on par- 
allel supercomputers, all have been realized. Efforts 
are currently underway to make GENRAL a parallel, 
multigrid FEM code. 

Implementation of iterative schemes like multigrid 
for solving the equations will make straightforward 
the application of accurate reflective and periodic 
boundary conditions. While possible with the Henyey 
technique, this process is much more difficult as it 
involves columns of matrix elements far from the di- 
agonal and special techniques for inversion. With it- 
erative methods, as with explicit codes, one can en- 
close the computational domain in a layer of "ghost" 
elements whose properties are determined at each it- 
eration by the interior solution. The ghost element 
approach will have the same order of accuracy as the 
interpolation scheme, unlike the current approach for 
the reflective boundary condition, which uses essen- 
tially a backward difference. 

Another possibly important issue is time evolution. 
All examples in this paper, even the four-dimensional 
ones, are time-independent and use a Cartesian met- 
ric. The inclusion of time dependence may be as sim- 
ple as employing a Minkowski metric and time deriva- 
tives of the variables, and the letting the finite element 
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machinery solve the problem. Often, however, the ad- 
dition of a new feature generates new numerical prob- 
lems which require modification of that machinery. 
Until more experience is obtained with time depen- 
dent problems, it is not clear whether the techniques 
discussed here are complete or whether they will need 
additional major development to handle evolutionary 
situations. 

Finally, many issues remain in the use of the finite 
element method for dynamical evolution problems. 
These are currently important topics in the engineer- 
ing field, but, because explicit finite difference codes 
do well for astrophysical problems in this area, devel- 
opment of these issues here will have lower priority. 
They include adaptive mesh refinement (for dynam- 
ical collapse situations), implementation of the gen- 
eral boundary conditions in equation (|^) (for magne- 
tohydrodynamics and solving Maxwell's or Einstein's 
equations), and proper upwinding schemes with be- 
havior comparable to the higher order Godonov schemes 
(for problems that develop shocks). 

5.4. A Note on Numerical Relativity with Fi- 
nite Element Analysis 

In numerical relativity it is customary to perform 
a "3-1-1 split" of the metric such that 



ds^ = (-a^ + (3i(3')dt^ + 2p,dx'dt +jijdx'dx^ 

(91) 

where 7ij is the 3-metric that raises or lowers indices 
on the shift 3- vector (3i, and a is the "lapse func- 
tion" , all of which are functions of position in space- 
time dArnowitt, Deser, fc Misner 196^ ; [York 19791 ). 



The 3-metric is specified on the initial hypersurface 
by solving the field constraints (initial value data), 
the lapse and shift are computed from four coordi- 
nate (or "gauge") conditions, and the Einstein field 
equations are used to evolve 7y to the next hyper- 
surface. The goal is to choose a gauge in which the 
hypersurfaces do not intersect a singularity before a 
significant amount of evolution occurs in some part of 
the mesh. The current method for singularity avoid- 
ance is to eliminate pathological parts o f spacetime 
from the mesh ("excise the black hole") ( Cook et al. 
19981) . 



Wh ile such an approach is also possib le with the 
FEM ( [Arnold, Mukherjee, fc Pouly 1998| , advancing 
time in a step-by-step fashion, the full covariant na- 
ture of the method and the lifting of the degeneracy 
between basic and mesh coordinate systems, allow ad- 



ditional approaches to be taken. In particular, it be- 
comes possible to extend the mesh fully in the time 
dimension, from initial to final hypersurface, choos- 
ing a relatively simple gauge for a and (3i. Then, 
adaptive gridding in all four dimensions can be used 
to keep the grid boundaries away from singularities 
and to further adjust the separation in time between 
spacelike hypersurfaces. Because of the ispoarametric 
transformation, the foliation no longer has to be along 
surfaces of constant time . The separation between 
adjacent surfaces can be non-uniform, the time coor- 
dinate can vary considerably over the hypersurfaces, 
and the final hypersurface even can end at different 
times. In effect, the adaptive gridding completes, in 
mesh coordinates, the job of slicing and singularity- 
avoiding that a poor gauge choice may fail to do. One 
advantage of this approach is that some or all of the 
field constraints can be applied on the final, instead 
of initial, hypersurface, turning an explicit hyperbolic 
problem into an implicit boundary value problem (like 
stellar structure) and possibly stabilizing the growth 
of errors. 

However, while such techniques probably can suc- 
ceed in keeping physical singularities at bay, it is 
doubtful that they can avoid coordinate singularities 
in general situations. (These arise in the most benign 
of curved surfaces — on the surface of the earth, for 
example.) Apart from knowing the geometry before 
hand and choosing the proper basic coordinate sys- 
tem, there are only a few ways to avoid these problems 
entirely. One is to embed the spacetime in a higher 
dimensional, flat Minkowskian space. In principle, 
as many dimensions as independent spacetime met- 
ric coefficients (ten) would be needed for the embed- 
ding, although it may be possible with fewer. From 
the ten hyperspace coordinates, and how they vary in 
the four-dimensional mesh, one then could derive the 
local metric of the spacetime and use it in the phys- 
ical equations. These ten equations for the metric in 
terms of the hyper-coordinates, plus the six Einstein 
equations and the four adaptive gridding equations, 
would be sufficient to determine the twenty indepen- 
dent gaf3 and hyper-coordinates at each node in the 
finite element mesh. While a fairly immense job for 
present-day computers, this prescription has the ad- 
vantage of being singularity-free in general situations. 

Another method of avoiding coordinate singulari- 
ties using the FEM is to dispense with global coor- 
dinates entirely, using only line segment lengths and 
deficit angles to describe the geometry and the Regge 
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calculus to describe the physics ( [Regge 1961 ; Hoist 
19£|8|). At present, this approach has been developed 
only for simplex-type elements and not hypercubes, 
so it is not a straightforward application of the code 
discussed herein. However, it may be useful to recast 
the Regge calculus for other element types. 

Finally, all methods that involve a full four-dimen- 
sional finite element spacetime are probably well be- 
yond the capabilities of present computer technology, 
even with the use of parallel multigrid techniques. 
Nevertheless, they appear to have such attractive fea- 
tures and elegance that it is important to begin to 
develop them. 
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A. Appendix A. Casting of the Differential 
Equations of Continuum Astrophysics into 
General Finite Element Form 

This appendix shows that virtually all the equa- 
tions of astrophysics of continuous media can be cast 
into the flux-conservative, finite element form (equa- 
tion |l]) and their boundary conditions into equation 
(||). That is, while possibly second order in spatial 
and time derivatives, they can be written as the four- 
divergence of a generalized stress tensor plus a gen- 
eralized body force vector, each of which are func- 
tions of no more than the first spacetime derivative 
(four-gradient) of the variables. Of course, it is al- 
ways possible to define additional variables (e.g., the 
24 connection coefficients) and turn the field equa- 
tions and conservation laws into first-order equations 
involving only the Fq term in equation (|^). The 



challenge, however, is to use only the original met- 
ric and field components as variables (avoiding addi- 
tional computational expense), and still maintain the 
flux-conservative form. Below is one solution to this 
problem. 

A.l. The Equations in Geometric Form 

The discussion here is concerned only with differen- 
tial equations. Local physics, such as the equations of 
state, opacity, emissivity, viscosity, etc., is not treated 
in detail. While having position and time dependence, 
these processes can be described with simple algebraic 
equations that do not affect the numerical method 
used. 

The differential equations are the deceptively sim- 
ple set of the Einstein equations for the gravitational 
field 

g = SttT (Al) 

(where Q and T are the symmetric Einstein cur- 
vature and stress-energy-momentum [SEM] tensors), 
and Maxwell's equations for the electromagnetic field 



V • :F = At: J 
V ■ M = 



(A2) 
(A3) 



where V is the covariant gradient operator, the an- 
tisymmetric Maxwell tensor Al = is the dual of 
the antisymmetric Faraday tensor and JT" is the 



four-current. With the symmetry, equations (Al) are 
10 in number, and (A2)-(A3) constitute 8, for a total 
of 18. However, because of identities satisfied by the 
Einstein and Faraday tensors, there are actually only 
12 independent equations (6 metric and 6 electromag- 
netic) but 16 unknowns at each point in space: the 
10 independent components of the metric g and the 
6 independent components (the electric and magnetic 
fields E and B) of the antisymmetric Faraday tensor. 
The remaining 4 metric unknowns are determined by 
the choice of a coordinate system or gauge. 

The standard method for generating a set of 12 
evolutionary equations is to project (A1)-(A3) into 
the hypersurface normal to a time-like vector (or 
world line) n with the projection tensor (equation 
For example, if — go^/ yZ—goo), then only the 
spatial part Sij will be non-zero, with i,j — 1,2,3. 
In general, however, n can be any time-like vector, 
so the equations will be left in general form. If the 



twelve spacelike components of equations (A1)-(A3) 
are satisfied throughout the four-dimensional space- 
time domain Q (with one factor of <S for each tensor 
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order) 

S-g-S = 8ttS-T-S (A4) 

S-{V-T) = inS-J (A5) 

S-{V-M) = (A6) 

then all that is necessary to satisfy the timelike com- 
ponents 

n-g ^ Sirn-T (A7) 
n-V-:F = Attu-J (A8) 
n-V-At = (A9) 

throughout all spacetime is to satisfy the latter equa- 
tions on one hypersurface only. Equations (A4)-(A£), 
therefore, are the 12 independent differential equa- 
tions to be solved for the six metric and six elec- 



tromagnetic field components, while equations ( A7)- 
( A9) are the constraints that need to be satisfied in or- 
der f or a solution to exist at all. (For reference, equa- 
tion (A5) is Ampere's law, ( [A^ ) Faraday's law, ( A7) 
contains the Hamiltonian and momentum constraints 
[by further contraction with n or «S, respectively], 
( A8) is Coulomb's law, and ( A9) is the solenoidal con- 
dition on the magnetic field.) 



Equations (A4), with (A7) as initial conditions, 
constitute the Cauchy problem of general relativ- 
ity. Equations ( [A^ ) and (A9) are the covariant 
form of the Evans-Hawley "constrained transport" 
method for enforcing the solenoidal constraint in non- 



Haw- 



relativistic magnetohydrodynamics (Evans 
ley 1988). Equations (A5) and (A8) represent con- 



strained transport in the presence of sources. When 



(A4)-(A6) are solved as Cauchy problems, equations 
(A7)-(A£) are applied on the initial hypersurface. 
However, as our approach here is to relax the system 
for a /our-dimensional spacetime, instead of evolving 
a three-dimensional surface forward in time, they can 
be applied on any spacelike hypersurface. 

In addition to the field equations, there are conser- 
vation laws that follow from identities satisfied by the 
fields. The Einstein curvature tensor is constructed 
in such a way that V • ^ = 0, so the conservation of 
energy and momentum 



V-T = 



(AlO) 



must also hold from equation (|A1|). Similarly, as ^ 
satisfies V- (V • ^) — 0, then the four-current must 
also be conserved 



The field equations then are "closed" by expressing 
the SEM tensor and four-current in terms of the state 
variables, and solving the conservation laws of energy 
and momentum for those variables. For most conceiv- 
able astrophysical situations — including those with 
niulti-fiuid dynamics, electromagnetic fields and cur- 
rents, radiation, viscosity, and nuclear reactions — 
expressions for T and JT" involve terms with, at most, 
first-order derivatives of the state variables with re- 
spect to space or time. This is true even in situa- 
tions near black hole horizons where particle interac- 
tion and fluid flow time scales are comparable, and 
equations like Ohm's law, for example, are no longer 
valid. 

A final group of differential equations comes from 
forming the zeroth, first, and second moments of the 
Boltzmann-Vlasov equation for each particle species 
(j) (photons, nuclei, etc.). These determine each 
species' individual number density, velocity (including 
the peculiar drift velocity q^-'^), and internal energy 
per particle. The kinetic equations give rise to famil- 
iar processes like nuclear burning, radiative transport, 
viscosity, and electrical conductivity. Nevertheless, 
they all have the same "conservative" form, with the 
divergence of a term that involves (at most) first-order 
spatial derivatives of the state variables. For example, 
the zeroth moment of these kinetic equations yields 



where A'^'-'^ s the particle number density and c^-" is 



V -J = 



(All) 



(A12) 

M) 

the net particle creation rate due to nuclear reactions. 
For many astrophysical situations, particle conserva- 
tion will be the only one of these equations needed, 
the drift velocity and particle energies being deter- 
mined by the diffusion or other approximations. 

A. 2. The Equations in Component Form 

With no loss of generality one can choose to write 
the differential equations in a coordinate frame. In 
this case, the connection coefficients are given by 

r"/57 = ^ g"^ + 9tJ.j,p - gpf-.tJ.) (ai3) 

and their trace reduces to 

F'^^^ = (lnV^)_^ (A14) 
where g is the determinant of the metric 

g = det\\g\\ (A15) 
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This yields simple expressions for the gradients and 
divergences in curved space, and the field equations 
(A4, A_5, and A6) and conservation laws ( A10| , All, 
and A12) become 



1 



9),P 



V-9 



(A16) 

(A17) 

(A18) 

(A19) 

(A20) 

(A21) 



which are all in the conservative finite element form 
(||). Use has been made of the expression for the 
Einstein tensor 



Gaf} 



n. 



1 



Q/3 



in terms of the Ricci tensor 
1 



'R-aP — 



-9 



where the Ricci scalar is the contraction TZ = TZ^ 
and the SEM scalar is the contraction T = T^, 



Note that, as equations ( A16D are symmetric in i and 
j, they are six in number, and (A17) and ( AlSD are 
three equations each. The term inside the divergence 
in equation (A16) is the trace-free curvature — the 



"flux" of the metric; the "force" terms involve deriva- 
tives of the projection tensor and a matter source gen- 
erating the curvature. 

Boundary conditions on the conservation laws are 
generally of the simple form (||). Furthermore, it is 
well known that the Einstein field constraints do not 
depend on second derivatives with respect to time, 
and the Maxwell constraints do not depend on time 
derivatives at all. So, it should be possible to write 
these in the form (^, with at most a divergence on 
the boundary and a gradient normal to it. One of the 
more compact versions of this is 



87r(n"n^r„^--r)]y^ = (A22) 

{5",K(r'',,rv-rvrV + ^^^0.0) + 

5",,^(n'^r^„0 - n^^T-^,)}^ = (A23) 



^ = (A24) 



-g = (A25) 
where n'^ is now identified as the boundary nor- 



mal. Equations (A22) and (A23) are, respectively, 
the Hamiltonian and momentum constraints. 



In addition to the Einstein field equations (A16) 
four more coordinate conditions, and their possible 
boundary conditions, are needed to completely spec- 
ify the metric. As these are rather arbitrary, these 
are left unspecified for now. (Of course, if they are 
to be differential equations, they also must be in the 
form (0).) 



B. Appendix B. Conditions on the Isopara- 
metric Transformation 

Not all £"q' specified by equation (0) are valid 
transformations in all situations. Sometimes they 
must satisfy certain conditions, which are discussed 
below. Two of these are conditions on the new (ele- 
ment mesh) metric formed by the transformation and 
the third is a condition on the element spacing itself. 
However, only the first condition always needs to be 
satisfied for a stable mesh. The other two are neces- 
sary only under certain circumstances. 

B.l. The Jacobi Condition 

The most important of these conditions is that the 
Jacobian determinant of the transformation be non- 
zero 

C = det||£|| ^ 

However, since £^17 — g' , and since g < 0, this also 
means that the determinant of the new metric must 
be negative 

g' = detllg'll < (Bl) 

The Jacobi condition ( p3l| ) ensures that the unit vec- 
tors in the element (primed) system, while not neces- 
sarily orthogonal, nevertheless form an independent 
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coordinate system that has an arrow of time. It is 
absolutely necessary that this condition be satisfied 
in order that the mesh be well behaved. 

B.2. The Local-Lorentz Condition 

A second possible condition is that the element co- 
ordinate system have a locally Lorentz character ev- 
erywhere — i.e., that be the mesh time coordinate. 
This ensures that, if the spatial portion of the element 
mesh moves with time, the mesh velocity always will 
be less than the speed of light. This is important, 
however, only if the mesh is used as a frame of refer- 
ence for measuring physical quantities. 

There are various ways of ensuring the Lorentz na- 
ture of the transformation C"' a'- The safest and sim- 
plest way is to ensure that each unit vector in the 
new space satisfies the proper timelike or spacelike 
constraint. With the constraints that denotes the 
time dimension and that ea> • epi = ga'f}', these local 
Lorentz conditions become 



5o'0' 
gi'i' 



= C 



{B2) 
(B3) 



Inequality ( B2 ) is equivalent to demanding that each 
of the element sides in the direction be spacelike. 

A less stringent, but still sufficient, condition on 
giH' could be derived by choosing a specific timelike 
vector, such as = Sq'^ /\/—go'0' (which still re- 
quires go'o' < 0), and then constructing from the cor- 
responding projection tensor a set of three indepen- 
dent vectors orthogonal to 

The condition that these vectors be spacelike {sji ■ 
Si' > 0) leads to a modified form for inequality ( |B^ ) 



> 



^(.9o'i'/\/-5o'0') 



In this less restrictive case, the e^/ unit vectors can 
be a bit timelike, but no more so than that given by 
the above inequality. 

B.3. The Courant Condition 

In standard nonrelativistic computational fluid dy- 
namics, in order that the (explicit) forward integra- 
tion in time be stable, the distance traversed in a 
single time step by sound waves or by the fluid itself 
(whichever is faster) must be substantially less than 



the mesh spacing. The ratio of these distances, called 
the Courant number C, is chosen to be ~ 0.1 — 0.4 or 
so depending on the stability of the numerical integra- 
tion scheme. This nonrelativistic Courant condition 
(which can be written as — At^/C^?;^^^. -|- Ax^ > 
for one-dimensional flow) easily generalizes in the 
four-dimensional general relativistic case to 

go.pAx'^Ax^ > 

where the vector Ax = {At/Cvmax, Ax, Ay, Az). It 
further generalizes in the case of a general curvilinear 
element mesh to 



gc.,p,AC' AC^' > 



(B4) 



in each element, with A^ = (A^*^ /Cw^^^., A^^ 



A^^ , A^^ ); A^" is the width of the element in each 
mesh dimension; and I'max is a three- velocity magni- 
tude equal to the maximum disturbance speed within 
the element. Because stable implicit techniques are 
used in time as well as space, a Courant number very 
close to unity probably can be tolerated. Therefore, 
in the case of relativistic flow, where v!^g_^ = 1, the 
Courant condition reduces to the requirement that 
the geodesies connecting opposing corners in each ele- 
ment must be spacelike (g^'/j/ A^" A^'' = As^ > 0). 



Generally, the Courant condition ( |B4D is much too 
restrictive and is routinely violated in slowly evolving 
or steady-state problems, where time steps are very 
long or even infinite. In implicit codes the condition 
needs to be satisfied only if one wishes to follow every 
short timescale transient phenomenon or wave. 
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FIGURE CAPTIONS 

Fig. 1. — Examples of shape functions in one dimen- 
sion for linear (top) and quadratic elements (bottom). 
The solid line shows the shape function for the node 
(at £,=£,() = 0)', the dashed line for the 1 node (at 
= 0.25); and so on. Note the element boundary 
nodes (large open circles) and interior nodes (smaller 
filled circles in the quadratic case). The derivatives 
of the shape functions are discontinuous at boundary 
nodes, although the functions themselves are contin- 
uous. Each function attains unit value at its corre- 
sponding node and exactly zero at all other nodes in 
the element. By definition, shape functions are also 
identically zero in elements not containing their cor- 
responding node. 

Fig. 2. — One-dimensional results for adaptive grid 
tests. Top left: solution of the Fermi-Dirac func- 
tion test with / = 50, using 8 linear, equally-spaced 
elements, and yielding a very large error (E^^ = 
0.27); top right: using 16 equally-spaced elements 
(E^2 = 0.027); middle left: 8 adaptive elements 
(E^^ = 0.031); middle right: 16 adaptive elements 
(E^2 ^ 0.0068); bottom left: same as middle right, 
but plotted against logw; bottom right: results when 
solving for w — log^Q w. The dotted line shows the 
exact solution. 

Fig. 3. — Two-dimensional results for adaptive grid 
tests, showing contour plots of the solution w. All 
use meshes of 16 x 16 linear elements with / = 50, 
a = 2, and 6 = c = in the Fermi-Dirac func- 
tion. In all cases, coordinates and differential equa- 
tions are expressed in x and y only, and although 
derivatives are calculated on the curvilinear grid, 
they are immediately transformed to the (x, y) sys- 
tem and used as such in the equations. Top left: 
uniformly-spaced Cartesian grid (E^^ = 0.081); top 
right: circular-polar grid (E^^ = 0.067); middle left: 
elliptical grid with the same axial ratio as the solu- 
tion (E^^ = 0.037); middle right: adaptive ellipti- 
cal grid allowing finer resolution of the Fermi surface 
(E^^ = 0.010); bottom left: same as middle right, but 
with logarithmic contours (note oscillations similar to 
those in the bottom left panel of Figure 2); bottom 
right: resulting solution when solving for w — log^o w- 



Fig. 4. — One-dimensional n — I spherical poly- 
tropic stars with 9 nodes. Left: standard Calerkin 
weighting of the FEM integrals, with linear elements 
(Ep^ = 1.4); middle: Petrov-Calerkin type 2 weight- 
ing and linear elements (E^^ — 0.048); right: same as 
middle, but with quadratic elements (E^^ = 0.0018). 
All models use adaptive gridding and logarithmic 
variables. Note the close nodal spacing near the stel- 
lar surface. 

Fig. 5. — Two-dimensional n = 1 spherical poly- 
tropic stars with 9^ nodes using the multipole bound- 
ary condition {w^'^). Top panels: mesh and pres- 
sure contours for (bi-) linear elements; bottom pan- 
els: mesh and pressure contours for (bi-) quadratic 
elements. Contours follow the precise interpolation 
within the respective elements. Note the adaptive 
gridding near the stellar surface, similar to that in 
Figure |^. 

Fig. 6. — Same as Figure |5[ but with 33^ nodes. Note 
the persistent errors in the pressure contours for linear 
elements which are absent in quadratic elements with 
far fewer nodes. 

Fig. 7. — Similar to Figure |6| but using the in- 
tegral boundary condition (wf) and solving a two- 
dimensional n = rotating polytropic star (Maclau- 
rin spheroid) with uP' j^-KGp = 0.224 — very close to 
the theoretical limit of 0.2246656. 

Fig. 8. — Comparison of rotating Maclaurin spheroid 
model stars, computed with different finite element 
meshes, with the analytic solution: Top panels: linear 
elements with 9^ and 17^ nodes, respectively; bottom 
panels: quadratic elements with 9^^ and 17^^ nodes. 
Although iJ^ is the primary model parameter, the fiat- 
tening ratio, total angular momentum, and uP' itself, 
with the normalizations in Tassoul (1978), are plotted 
as a function of the derived parameter r for compar- 
ison with the analytic theory. 

Fig. 9. — Absolute value of the fractional errors for 
the Maclaurin spheroids in Figure ^ as a function of 
LiP . Solid lines show models with linear elements, 
dashed lines quadratic elements; curves without sym- 
bols are for 9^ node models, those with are for 17^ 
models. 
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